#####################################################################################
## GATT/WTO and Trade application
## Code for replicating the three models reported in the article
############################################
## Initial Analysis of WTO and Free Trade
## Generate Figures A13 and A14(b)
###############################################
rm(list = ls())


sink("log/log_WTOReplicationInitial.txt") # start log
begin.time<-Sys.time()


## compile functions 
library(bpNet)
#library(mcmcr)
library(abind)
library(ggplot2)
library(ggstance)
library(coda)
library(igraph)
library(coda)
library(stargazer)

aa <- -0.6  ##??
bb <- 0.1    ## ??
set.seed(123456789)

## ------------------------- ##

## setwd('~/Desktop/PSRMreplicationLP/data/EmpData')
source('code/net_plot.R', chdir = TRUE)    ## load code to make figures

## Load data
traderep <- read.csv("data/EmpData/tradewlag.csv")
## set the missing data as NA rather than -999
traderep$l1yrsoffc[traderep$l1yrsoffc==-999] <- 0
rhoZ1 <- read.csv("data/EmpData/WTOdensity.csv")[,-1]
load("data/EmpData/WmWTOlist.RData") # WTO network

attach(traderep)

 Y <- newtar[year>1970] # because the WTO network starts from 1971 in the replication data
 n <- length(Y)
 X1 <- cbind(l1polity2, l1log.pop, l1log.gdppc,l1bpc,l1ecris,l1yrsoffc,
              l1gatt.wto, l1usheg, l1avetar, Eight, Wmember, wtrade, lagnewtar)[year>1970,]

 Uid <- cowcode[year>1970]
 Tid <- traderep$year[year>1970]
dataWTO <- data.frame(Uid, Tid, Y, X1)
## produce Table A3 (see file 5_WTOreplication.R)
## stargazer(dataWTO[,-c(1:2)])
 

 ## ------------- gen W matrix ------------ ## 
 W <- WmWTOlist[-1]
 W2 <- array(as.numeric(unlist(W)), dim=c(85, 85, 38))

 ## Get dimension of TSCS 
 N <- dim(W2)[1]
  TT <- dim(W2)[3]

## correct the W matrix (## The matrix has the diagnal nonzero)
## and standardize the matrix
 W3 <- W2
 for (i in 1:TT){
    w <-  W2[ , , i]
    diag(w) <- 0
    ww <- (w/rowSums(w))
    ww[ww=="NaN"] <- 0
    W3[ , , i] <- ww
    }

## combine the network measure and group-level intercept
rhoZ <- cbind(1, rhoZ1)

#int <- c(1, 1, rep(0, 36)) ## two 1's because we include y lag and the first period will be taken out
#rhoZ2 <- cbind(int, rhoZ1)

zero <- cbind(rep(0, 38))

## number of total iterations and number of burnin
mcmc <- 55000
burnin <- 5000

## Factor model
outMF <- bpNet(data = dataWTO,   
	              W = W3,         
	              index = c("Uid", "Tid"), 
	              Yname = "Y",
	              Xname = c("l1polity2", "l1log.pop", "l1log.gdppc","l1bpc","l1ecris","l1yrsoffc",  "l1gatt.wto", "l1usheg", "l1avetar", "Eight",  "wtrade"), 
	              Zname = NULL,  
	              Aname = NULL,   
                      Contextual = c("l1bpc","l1ecris"),  
	              Contextual.effect = "both", 
	              lagY = TRUE,                
	              force = "both",             
	              r = 30,          
	              flasso = 1,             
	              rhoZ = as.matrix(rhoZ),      
	              rho_spec = "multilevel",               
                      niter = mcmc,
              burn = burnin)



############### output and graph-making
## Figure A11: factor selection
p1 <- omegaPlot(outMF, oxlim = c(-5, 5), plotlength = 5, labelname = "Factor")
ggsave('plots/WTO/factors.pdf', p1, width = 18, height = 12)

## Figure A10(b): factor selection
rho <- rep(0, 37)
TT <- 37
modelRho <- outMF$Rho
		## ----------------------- 1 varying rho_t
datap <- cbind.data.frame(1972:2008, apply(modelRho, 1, mean),
		          apply(modelRho, 1, quantile, 0.025),
		          apply(modelRho, 1, quantile, 0.975),
		          c(rho))
names(datap) <- c("time", "mean", "cil", "ciu", "true")
datap$type <- rep("m", TT)
p <- ggplot(datap) + xlab("Time") + ylab("Rho")
	p <- p + geom_line(aes(time, mean, colour = type, size = type, linetype = type))
	p <- p + geom_line(aes(time, true), size = 0.5, colour = "red", linetype = "dashed")
	p <- p + geom_ribbon(aes(x = time, ymin = cil, 
	                                   ymax = ciu), fill = "#000000", alpha = 0.2)

	set.limits <- c("t", "m", "c")
	set.labels <- c("0", "Posterior Mean", "95% CI")
	set.colors <- c("red", "grey50", "#00000080")
	set.size <- c(1, 1, 3)
	set.type <- c("dashed", "dashed", "solid")

	p <- p + scale_colour_manual(limits = set.limits,
	                             labels = set.labels,
	                             values =set.colors) +
	         scale_size_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.size) +
	         scale_linetype_manual(limits = set.limits,
	                           labels = set.labels,
	                           values = set.type) +
	         guides(colour = guide_legend(title=NULL, nrow=1),
	                size = guide_legend(title=NULL, nrow=1),
	                linetype = guide_legend(title=NULL, nrow=1))
	## p <- p + theme_bw(base_size = 15)
	p <- p + theme(
	               legend.text = element_text(margin = margin(r = 10, unit = "pt"), size = 12),
	               legend.position = "bottom",
	               axis.title = element_text(size=12),
	               axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0)),
	               axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
	               axis.text = element_text(color="black", size=8),
	               axis.text.x = element_text(size = 8),
	               axis.text.y = element_text(size = 8),
	               plot.title = element_text(size = 15,
	                                         hjust = 0.5,
	                                         face="bold",
	                                         margin = margin(10, 0, 10, 0)))
p <- p + scale_y_continuous(limits = c(-1.2, 0.05))
a <- 16
p <- p +theme(axis.text.x = element_text( hjust = 1, size=a))
p <- p +theme(axis.text.y = element_text( hjust = 1, size=a))
ggsave('plots/WTO/RhoMF230.pdf', p, width = 9, height = 6)

print(Sys.time()-begin.time) 
sink() # close log  



